Numerical Methods for Integration

Quadratures

Introduction

A quadrature is just a name we use to describe methods of numerical integration

In this chapter we will look at some simple methods such as the trapezium rule, which you will already be familiar with and then we will consider some more sophisticated developments, such as Simpson's and Boole's rule, which interpolate function points with polynomials of increasing degree

Finally we will look at Gaussian quadratures where the assumption of evenly spaced abscissa points is dropped and a more general form of quadrature is developed which can produce remarkably accurate results with relatively little calculation

Very basic quadrature methods

The most basic form of numerical integration is called the left point rule: In this case we simply divide the function into intervals of width $\Delta x$ and then add up the area of the rectangles width: $\Delta x$ and height: $f(x_i)$, where $x_i$ is the value of $x$ at the left of the interval

This gives: $\int_{a}^{b}f(x) dx \approx \Delta x \times \left(f(x_0)+f(x_1)+...+f(x_{n-1})\right)$

where $a=x_0$, $b=x_n$ and the $x_i$s are evenly spaced across the interval.

The right point rule and the mid-point rule follow by analogy

Clearly these methods are very primitive and it is easy to improve on them

The Trapezium Rule

The trapezium rule represents the first level of sophistication in the numerical calculation of integrals

consider the problem of how to calculate the area under the curve: $y=e^{-x} \times sin x$ between $0$ and $\pi$

The following graph shows how we could use a series of trapeziums to estimate this integral

We can easily see that the area under the trapeziums reduces to:

$A=\frac{\pi}{8} \times \left(f(0) + 2f(\frac{\pi}{4}) + 2f(\frac{\pi}{2}) + 2f(\frac{3\pi}{4}) + f(\pi) \right)$

So more generally the trapezium rule is given by:

$\int_{a}^{b} f(x) dx \approx \frac{\Delta x}{2} \times \left(f(x_0)+2f(x_1)+2f(x_2)+...+2f(x_{n-1})+f(x_{n})\right)$

Simpson's Rule

Suppose instead of fitting straight line segments to our function we attempt to find a better fit by fitting quadratic segments

Look at the same function below (grey) with two black quadratics fitted to the first three and last three calculated function points

This leads us to Simpson's Rule which is stated generally as:

$\int_{a}^{b} f(x) dx \approx \frac{\Delta x}{3} \times \left(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+...+2f(x_{n-2})+4f(x_{n-1})+f(x_{n})\right)$

Simpson's three eighths Rule

The next level is to fit cubic segments

Look at the same function again (grey) with the red cubic fitted to four calculated function points

The principles are exactly the same but the algebra is a bit more complicated, so Simpson's 3/8ths Rule over a $3\Delta x$ interval can be expressed as:

$\int_{a}^{b} f(x) dx \approx \frac{3\Delta x}{8} \times \left(f(x_0)+3f(x_1)+3f(x_2)+f(x_3)\right)$

Boole's Rule

The next level is to fit quartic segments

This time the interpolating function is a quartic

The principles are again the same, so Boole's Rule over $4\Delta x$ interval can be expressed as:

$\int_{a}^{b} f(x) dx \approx \frac{\Delta x}{45} \times \left(14f(x_0)+64f(x_1)+24f(x_2)+64f(x_3)+14f(x_4)\right)$

Programming Extensions

There are a number of things you should now do to tidy up your routine

• Ensure the routine checks for an appropriate number of intervals for each integration method
• You could splice different methods where $n$ does not fit
• You could include an option to leave out $n$ and let the routine try different values until some convergence criterion is reached

Gaussian Quadratures

So far all of our methods have assumed that the abscissa points have been equally spaced along the $x$-axis

We could have changed this if we wanted but there seemed little point and it would have made the algebra more complicated

The question therefore arises: Can we improve the speed or accuracy of our routine by using a different set of abscissa points?

The answer is YES and the set of methods we can use are referred to as Gaussian quadrature

General Considerations

In general all quadrature methods are based on the formula:

$\displaystyle\int_{a}^{b}f(x)dx\approx \displaystyle\sum_{i=1}^{n}\omega_i f(x_i)$

Where $\omega_i$ are a set of weights and $x_i$ are a set of abscissa points

As we are now free to choose $2n$ values we should be able to approximate our integral much more accurately - in fact if we motivate our selection of $\omega_i$ and $x_i$ by fitting polynomials then this method will calculate the integral of an order $2n-1$ polynomial exactly

2 Point Gaussian Quadratures

W.l.o.g we work on the interval $[-1, 1]$

We wish to find $\omega_1$, $\omega_2$, $x_1$ and $x_2$ such that

$\displaystyle\int_{-1}^{1}f(x)dx \equiv \displaystyle\sum_{i=1}^{2}\omega_i f(x_i)$, for a cubic $f(x)$, i.e.:

$\displaystyle\int_{-1}^{1} a_0+a_1x+a_2x^2+x_3x^3dx \equiv \omega_1 f(x_1) + \omega_2 f(x_2)$, i.e. this relation must hold for all choices of $a_i$

We proceed with basic calculus

$\left[ a_0x+\frac{a_1x^2}{2}+\frac{a_2x^3}{3}+\frac{x_3x^4}{4} \right]_{-1}^{1} \equiv \omega_1 f(x_1) + \omega_2 f(x_2)$

$a_0(1-(-1)) + a_1\frac{1^2-(-1)^2}{2} + a_2\frac{1^3-(-1)^3}{3} + a_3\frac{1^4-(-1)^4}{4} \equiv \omega_1 f(x_1) + \omega_2 f(x_2)$

$2a_0 + \frac{2}{3}a_2 \equiv \omega_1 \left(a_0+a_1x_1+a_2x_1^2+a_3x_1^3\right)+ \omega_2\left(a_0+a_1x_2+a_2x_2^2+a_3x_2^3\right)$

$2a_0 + \frac{2}{3}a_2 \equiv a_0(\omega_1 + \omega_2) + a_1(\omega_1 x_1 + \omega_2 x_2) + a_2(\omega_1 x_1^2 + \omega_2 x_2^2) + a_3(\omega_1 x_1^3 + \omega_2 x_2^3)$

As this is an identity we can match the co-efficients of the $a_i$

$a_0$: $2 = \omega_1 + \omega_2$

$a_1$: $0 = \omega_1 x_1 + \omega_2 x_2$

$a_2$: $\frac{2}{3} = \omega_1 x_1^2 + \omega_2 x_2^2$

$a_3$: $0 = \omega_1 x_1^3 + \omega_2 x_2^3$

The best way to solve these equations is guess and check as follows:

as symmetry seems likely assume $\omega_1 = \omega_2$ therefore $\omega_1 = \omega_2 = 1$

$a_1$ gives $x_1 = -x_2$

$a_2$ gives $\frac{2}{3} = \omega_1 x_1^2 + \omega_2 x_2^2 = 2 x_1^2$

so $3x_1^2-1=0$ - this is the Legendre Polynomial order 2

w.l.o.g. $x_1=-\sqrt{\frac{1}{3}}$ and $x_2=\sqrt{\frac{1}{3}}$

The last thing we need to do to make our method useful is to drop the requirement to work on the $[-1, 1]$ interval

There are two additional things we need to consider to to this:

• the $\omega_i$ need to be transposed from the [-1,1] interval to the [5,10] interval
• The area calculated at the end needs to be ratioed up from a width of 2 (1 - (-1)) to a width of 5 (10-5)

n Point Gaussian Quadratures

We could repeat the above exercise for 3 and then 4,5 ,6 etc points. The algebra does get complicated, but fortunately we have the following theorem:

We do not attempt to prove this result here - but the above analysis is the proof of this result for 2 point quadrature. You can see how the Legendre Polynomial arises in the above analysis

The Legendre Polynomials are defined by: $P_n(x)=\frac{1}{2^n} \displaystyle\sum_{k=0}^{n} \left({n \choose k}^2 \times (x-1)^{n-k} \times (x+1)^k\right)$

Formula 98 files are Legendre Polynomials and Legendre Polynomials wi

Calculating $P'_n(x)$ is a simple application of the product rule

So we can see that the first few of the Legendre Polynomials are:

• $n=0: P_0(x)=0$
• $n=1: P_1(x)=x$
• $n=2: P_2(x)=\frac{1}{2} \left(3x^2-1 \right)$
• $n=3: P_3(x)=\frac{1}{2} \left(5x^3-3x \right)$
• $n=4: P_4(x)=\frac{1}{8} \left(35x^4-30x^2 + 3 \right)$

How might we use the formula: $x=cos \left( \frac{(i-0.25) \pi}{n+0.5} \right)$

To produce initial estimates of the roots so we can then use Newton Rhapson to produce our exact roots

Fourier Analysis

Fourier analysis allows us to convert waveforms into spectral densities and back again

You can soon see that as $n \rightarrow \infty$ this function becomes a sawtooth wave

This raises the question: If we started with the sawtooth wave wave could we calculate the co-efficients

This is where fourier analysis comes in

For our purposes we will skip Fourier series and go onto the continuous Fourier transform

For our purposes (later in the course when we look at stochastic volatility) we use Fourier inversion to get from the characteristic function of a probability distribution back to the probability distribution

The above formulas look rather daunting but we start by breaking them down with Euler's formula:

$e^{ix}=cos x + i sin x$

Investigating whether or not the value calculated above actually makes any sense will form part of your written coursework. For the calculation part of the coursework it is sufficient to be able to just calculate what the actual number is.